#delimit;
clear;

use 20161021_Working_hh.dta, clear;
keep if year>=2011 & year<=2015;
drop if statefip==11;
rename statefip FIPS;
merge m:1 FIPS using "FIPS.dta", nogen;


replace Y=(Y/1000)*CPIU2016;
replace totalTransfers=(totalTransfers/1000)*CPIU2016;

save 20161021_Working_hh_1014.dta, replace;

clear;
save 20161021_parametersByState.dta, replace emptyok;

local list "AK AL AR AZ CA CO CT DE FL GA HI IA ID IL IN KS KY LA MA MD ME MI MN MO MS
MT NC ND NE NH NJ NM NV NY OH OK OR PA RI SC SD TN TX UT VA VT WA WI WV WY" ;


foreach i of local list  {;
matrix drop _all;
use 20161021_Working_hh_1014.dta, clear;
  nl (totalTransfers = {a} + {b1}*exp({b2}*Y)) if state_al=="`i'", initial(a 1 b1 2 b2 -.1) vce(gnr);
  matrix B=e(b);
  matrix V=e(V);
  matrix v=(V[1,1], V[2,1], V[2,2], V[3,1] ,V[3,2], V[3,3]);
  
	
  clear;
  svmat B;
  svmat v;
  gen state_alpha="`i'";
  append using 20161021_parametersByState.dta;
  save 20161021_parametersByState.dta, replace;

};


use FIPS.dta, clear;
sort state_alpha;

merge 1:1 state_alpha using 20161021_parametersByState.dta ;
drop if _merge~=3;
drop _merge;

rename B1 A;
rename B2 B1;
rename B3 B2; 
rename v1 V_A;
rename v2 C_AB1;
rename v3 V_B1;
rename v4 C_AB2;
rename v5 C_B1B2;
rename v6 V_B2; 


sort state;
save 20161021_parametersByState, replace;



*>>>>A  Constructing Matrices;

*/;

use 20161021_parametersByState, clear;

matrix drop _all;

mkmat A B1 B2, matrix(BETA) rown(state_alpha);



mkmat V_A C_AB1 V_B1 C_AB2 C_B1B2 V_B2  , matrix(COV) rown(state_alpha);



local list "AK AL AR AZ CA CO CT DE FL GA HI IA ID IL IN KS KY LA MA MD ME MI MN MO MS
MT NC ND NE NH NJ NM NV NY OH OK OR PA RI SC SD TN TX UT VA VT WA WI WV WY" ;


foreach i of local list  {;

    matrix BETA_`i'= BETA["`i'", 1..3];
    matrix COV_`i'=COV["`i'", 1..6];


        
};


*>>>>B  Generating Simulated Parameters;

drop _all;
save "20161021_Simulations.dta", replace emptyok;


local list "AK AL AR AZ CA CO CT DE FL GA HI IA ID IL IN KS KY LA MA MD ME MI MN MO MS
MT NC ND NE NH NJ NM NV NY OH OK OR PA RI SC SD TN TX UT VA VT WA WI WV WY" ;


foreach i of local list  {;


    drawnorm ALPHA BETA_1 BETA_2, n(1000) means(BETA_`i') cov(COV_`i') cstorage(lower) seed(00001) clear;
    generate state_alpha="`i'";
    
	append using "20161021_Simulations.dta";
    save "20161021_Simulations.dta", replace;

 
    
};

*>>>>C  Generating Distribution of TAU and R;


use 20161021_Simulations.dta, clear;
sort state;

merge m:1 state_alpha using 20161021_parametersByState; 


local list "AK AL AR AZ CA CO CT DE FL GA HI IA ID IL IN KS KY LA MA MD ME MI MN MO MS
MT NC ND NE NH NJ NM NV NY OH OK OR PA RI SC SD TN TX UT VA VT WA WI WV WY" ;

gen PSI_1=.5*19.790;
gen PSI_2=1*19.790;
gen PSI_3=1.5*19.790;
gen PSI_4=2*19.790;

forvalues p=1/4{;

generate TAU_`p'=.7*PSI_`p';

foreach i of local list  {;

local j=1;

while `j'  < 100  {;
    generate TAU_`p'_`i'_`j'= TAU_`p' +`j'*.1 if (state_al=="`i'") ;
	generate TEST_`p'_`i'_`j'=(PSI_`p'-TAU_`p'_`i'_`j') - (ALPHA+BETA_1*exp(BETA_2*TAU_`p'_`i'_`j')) if (state_al=="`i'") ;
	replace TAU_`p'=TAU_`p'_`i'_`j' if TEST_`p'_`i'_`j'<.1 & TEST_`p'_`i'_`j'>0;
	
	local j=`j'+1;
	};



local k=1;

while `k'  < 100  {;
    replace TAU_`p'_`i'_`k'= TAU_`p' - `k'*.01 if (state_al=="`i'") ;
	replace TEST_`p'_`i'_`k'=(PSI_`p'-TAU_`p'_`i'_`k') - (ALPHA+BETA_1*exp(BETA_2*TAU_`p'_`i'_`k')) if (state_al=="`i'") ;
	replace TAU_`p'=TAU_`p'_`i'_`k' if TEST_`p'_`i'_`k'<.01 & TEST_`p'_`i'_`k'>0;
	
	local k=`k'+1;
	};


local l=1;

while `l'  < 100  {;
    replace TAU_`p'_`i'_`l'= TAU_`p' - `l'*.001 if (state_al=="`i'") ;
	replace TEST_`p'_`i'_`l'=(PSI_`p'-TAU_`p'_`i'_`l') - (ALPHA+BETA_1*exp(BETA_2*TAU_`p'_`i'_`l')) if (state_al=="`i'") ;
	replace TAU_`p'=TAU_`p'_`i'_`l' if TEST_`p'_`i'_`l'<.001 & TEST_`p'_`i'_`l'>0;
	
	local l=`l'+1;
	};



drop TEST* TAU_`p'_*;

};

generate R_`p'= ((A*TAU_`p' + (B1/B2)*exp(B2*TAU_`p') - (B1/B2))  + (PSI_`p'^2 - PSI_`p'^2/2) - (PSI_`p'*TAU_`p' - TAU_`p'^2/2)) /(PSI_`p'^2 - PSI_`p'^2/2);


};
                                                                                  

/* NB Corrected equation */;     


#delimit;
collapse A B1 B2  V_A V_B1 V_B2 TAU* PSI* R* C* (sd) SD_R_1=R_1 SD_R_2=R_2 SD_R_3=R_3 SD_R_4=R_4 SD_TAU_2=TAU_2, by(FIPS); 

merge 1:1  FIPS using "FIPS.dta", nogen;


save "20161021_State_Parameters_w_Var.dta", replace;








